Remigio-Capstone-InSitu

library(dplyr)
library(raster)
library(rhdf5)
library(rgdal)
library(neonAOP)
library(reshape2)
library(ggplot2)
library(MASS)

options(stringsAsFactors = FALSE)

Import D17-SOAP Field Data

## read in file
insitu.SOAP<-read.csv("/Users/R-Files/Documents/data/NEONDI-2016/NEONdata/D17-California/SOAP/2013/insitu/veg-structure/D17_2013_SOAP_vegStr.csv")

## visualize structure
str(insitu.SOAP)

## 'data.frame':    1231 obs. of  30 variables:
##  $ siteid               : chr  "SOAP" "SOAP" "SOAP" "SOAP" ...
##  $ sitename             : chr  "Soaproot Saddle" "Soaproot Saddle" "Soaproot Saddle" "Soaproot Saddle" ...
##  $ plotid               : chr  "SOAP43" "SOAP43" "SOAP43" "SOAP43" ...
##  $ easting              : num  297065 297070 297056 297067 297057 ...
##  $ northing             : num  4100726 4100722 4100722 4100719 4100729 ...
##  $ taxonid              : chr  "ABCO" "ABCO" "ABCO" "ABCO" ...
##  $ scientificname       : chr  "Abies concolor" "Abies concolor" "Abies concolor" "Abies concolor" ...
##  $ indvidual_id         : int  622 632 567 606 639 566 527 638 516 555 ...
##  $ pointid              : chr  "center" "NE" "center" "center" ...
##  $ individualdistance   : num  5.4 5.4 4.1 8.5 5.7 4.1 8.3 6.9 6.4 6 ...
##  $ individualazimuth    : num  67.8 188.7 250.8 122.4 332.6 ...
##  $ dbh                  : num  67.6 16.1 5.9 0.9 4 5.1 1.1 3.1 1 8.9 ...
##  $ dbhheight            : num  130 130 130 10 130 130 130 130 10 130 ...
##  $ basalcanopydiam      : num  8.2 0 0 0 0 0 0 0 0 0 ...
##  $ basalcanopydiam_90deg: num  7.9 0 0 0 0 0 0 0 0 0 ...
##  $ maxcanopydiam        : num  0 3.3 2.9 0.4 2 2.4 1.3 2 0.5 3 ...
##  $ canopydiam_90deg     : num  0 3.2 2.8 0.3 1.5 1.9 1.4 1.9 0.5 2.9 ...
##  $ stemheight           : num  34.5 13 3.8 0.6 3.4 3 1.4 3.1 0.5 5 ...
##  $ stemremarks          : chr  "" "" "" "" ...
##  $ stemstatus           : chr  "" "" "" "" ...
##  $ canopyform           : chr  "" "" "" "" ...
##  $ livingcanopy         : int  100 100 100 100 100 100 100 100 100 100 ...
##  $ inplotcanopy         : int  100 100 100 100 100 100 100 100 100 100 ...
##  $ materialsampleid     : chr  "" "" "" "" ...
##  $ dbhqf                : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ stemmapqf            : int  0 1 0 0 0 0 0 0 0 0 ...
##  $ plant_group          : logi  NA NA NA NA NA NA ...
##  $ common_name          : logi  NA NA NA NA NA NA ...
##  $ aop_plot             : logi  NA NA NA NA NA NA ...
##  $ unique_id            : logi  NA NA NA NA NA NA ...

#summarize
summary(insitu.SOAP)

##     siteid            sitename            plotid             easting      
##  Length:1231        Length:1231        Length:1231        Min.   :297034  
##  Class :character   Class :character   Class :character   1st Qu.:297064  
##  Mode  :character   Mode  :character   Mode  :character   Median :298705  
##                                                           Mean   :298489  
##                                                           3rd Qu.:299809  
##                                                           Max.   :300153  
##                                                                           
##     northing         taxonid          scientificname      indvidual_id   
##  Min.   :4100083   Length:1231        Length:1231        Min.   :   1.0  
##  1st Qu.:4100556   Class :character   Class :character   1st Qu.: 304.5  
##  Median :4100844   Mode  :character   Mode  :character   Median : 612.0  
##  Mean   :4100817                                         Mean   : 632.2  
##  3rd Qu.:4101027                                         3rd Qu.: 920.5  
##  Max.   :4101489                                         Max.   :1341.0  
##                                                                          
##    pointid          individualdistance individualazimuth      dbh        
##  Length:1231        Min.   : 0.000     Min.   :  0.00    Min.   :  0.00  
##  Class :character   1st Qu.: 4.700     1st Qu.: 66.95    1st Qu.:  1.20  
##  Mode  :character   Median : 6.900     Median :160.80    Median :  3.65  
##                     Mean   : 6.903     Mean   :173.57    Mean   : 10.10  
##                     3rd Qu.: 9.100     3rd Qu.:286.30    3rd Qu.: 13.68  
##                     Max.   :64.000     Max.   :359.80    Max.   :159.80  
##                                                          NA's   :249     
##    dbhheight     basalcanopydiam   basalcanopydiam_90deg maxcanopydiam   
##  Min.   :  0.0   Min.   :  0.000   Min.   : 0.000        Min.   : 0.000  
##  1st Qu.: 10.0   1st Qu.:  0.000   1st Qu.: 0.000        1st Qu.: 0.700  
##  Median :130.0   Median :  0.000   Median : 0.000        Median : 1.600  
##  Mean   : 86.3   Mean   :  7.556   Mean   : 5.322        Mean   : 2.355  
##  3rd Qu.:130.0   3rd Qu.:  4.000   3rd Qu.: 2.000        3rd Qu.: 3.100  
##  Max.   :140.0   Max.   :220.000   Max.   :93.000        Max.   :80.000  
##                  NA's   :1                               NA's   :1       
##  canopydiam_90deg   stemheight      stemremarks         stemstatus       
##  Min.   : 0.000   Min.   : -1.600   Length:1231        Length:1231       
##  1st Qu.: 0.500   1st Qu.:  1.000   Class :character   Class :character  
##  Median : 1.300   Median :  1.800   Mode  :character   Mode  :character  
##  Mean   : 1.863   Mean   :  4.786                                        
##  3rd Qu.: 2.400   3rd Qu.:  4.450                                        
##  Max.   :70.000   Max.   :134.000                                        
##  NA's   :1                                                               
##   canopyform         livingcanopy     inplotcanopy    materialsampleid  
##  Length:1231        Min.   :  0.00   Min.   : 40.00   Length:1231       
##  Class :character   1st Qu.: 60.00   1st Qu.:100.00   Class :character  
##  Mode  :character   Median :100.00   Median :100.00   Mode  :character  
##                     Mean   : 74.25   Mean   : 99.78                     
##                     3rd Qu.:100.00   3rd Qu.:100.00                     
##                     Max.   :100.00   Max.   :100.00                     
##                     NA's   :85       NA's   :84                         
##      dbhqf     stemmapqf       plant_group    common_name   
##  Min.   :0   Min.   :0.00000   Mode:logical   Mode:logical  
##  1st Qu.:0   1st Qu.:0.00000   NA's:1231      NA's:1231     
##  Median :0   Median :0.00000                                
##  Mean   :0   Mean   :0.09748                                
##  3rd Qu.:0   3rd Qu.:0.00000                                
##  Max.   :0   Max.   :1.00000                                
##  NA's   :1                                                  
##  aop_plot       unique_id     
##  Mode:logical   Mode:logical  
##  NA's:1231      NA's:1231     
##                               
##                               
##                               
##                               
## 

Exploratory Analysis: Subset variables of interest and identify data type (Cat & Cont)

#  selected potential predictors of interest at SOAP 2013

# plotid, taxonid, scientificname, dbh, dbhheight, basalCanopyDiam, Maximum canopy diameter, stemheight
# canopyform, livingcanopy,  inplotcanopy

# subset dataframe-- make it not dirty
insitu.SOAP<-subset(insitu.SOAP, select=c(plotid, easting, northing, taxonid, scientificname, dbh, dbhheight, basalcanopydiam, maxcanopydiam, stemheight, canopyform, livingcanopy,  inplotcanopy))

## subset categorical variables
SOAP.cat<-subset(insitu.SOAP, select=c(plotid, taxonid, scientificname, canopyform))

## subset continuous variables
SOAP.cont<- subset(insitu.SOAP, select=c(dbh, dbhheight, basalcanopydiam, maxcanopydiam, stemheight, livingcanopy, inplotcanopy))

Exploratory Analysis for SOAP- Frequency Categorical Statistics

## frequency summaries by:

# plotid
SOAP.cat %>%
  group_by(plotid) %>%
  summarise (n = n()) %>%
  mutate(rel.freq = paste0(round(100 * n/sum(n), 1), "%"))

## Source: local data frame [18 x 3]
## 
##      plotid     n rel.freq
##       (chr) (int)    (chr)
## 1  SOAP1343    24     1.9%
## 2   SOAP139   122     9.9%
## 3   SOAP143   114     9.3%
## 4  SOAP1515    26     2.1%
## 5  SOAP1563    14     1.1%
## 6  SOAP1611    10     0.8%
## 7  SOAP1695   106     8.6%
## 8   SOAP187    95     7.7%
## 9   SOAP223    43     3.5%
## 10  SOAP255   104     8.4%
## 11  SOAP283    14     1.1%
## 12  SOAP299    26     2.1%
## 13  SOAP331   132    10.7%
## 14   SOAP43   100     8.1%
## 15  SOAP555     6     0.5%
## 16   SOAP63   122     9.9%
## 17   SOAP95   106     8.6%
## 18  SOAP991    67     5.4%

# dominant(?) plant type

SOAP.cat %>%
  group_by(scientificname) %>%
  summarise (n = n()) %>%
  mutate(rel.freq = paste0(round(100 * n/sum(n), 1), "%"))

## Source: local data frame [16 x 3]
## 
##                      scientificname     n rel.freq
##                               (chr) (int)    (chr)
## 1                    Abies concolor    11     0.9%
## 2             Amelanchier utahensis     2     0.2%
## 3            Arctostaphylos viscida   179    14.5%
## 4              Calocedrus decurrens   435    35.3%
## 5                Ceanothus cuneatus    49       4%
## 6            Ceanothus integerrimus    42     3.4%
## 7  Cercocarpus montanus var. glaber    40     3.2%
## 8                   Corylus cornuta     3     0.2%
## 9                 Pinus lambertiana    16     1.3%
## 10                  Pinus ponderosa   213    17.3%
## 11                Prunus emarginata     1     0.1%
## 12              Quercus chrysolepis   141    11.5%
## 13                 Quercus garryana     3     0.2%
## 14                Quercus kelloggii    65     5.3%
## 15               Rhamnus ilicifolia    12       1%
## 16                    Ribes roezlii    19     1.5%

# canopyform

SOAP.cat %>%
  group_by(canopyform) %>%
  summarise (n = n()) %>%
  mutate(rel.freq = paste0(round(100 * n/sum(n), 1), "%"))

## Source: local data frame [6 x 3]
## 
##   canopyform     n rel.freq
##        (chr) (int)    (chr)
## 1              844    68.6%
## 2       Cone   155    12.6%
## 3   Cylinder    14     1.1%
## 4 Hemisphere    70     5.7%
## 5     Sphere    64     5.2%
## 6         NA    84     6.8%

Exploratory Analysis for SOAP- Summary Continuous Statistics

## Continuous variable summaries- Min/Max, Quants
summary(SOAP.cont)

##       dbh           dbhheight     basalcanopydiam   maxcanopydiam   
##  Min.   :  0.00   Min.   :  0.0   Min.   :  0.000   Min.   : 0.000  
##  1st Qu.:  1.20   1st Qu.: 10.0   1st Qu.:  0.000   1st Qu.: 0.700  
##  Median :  3.65   Median :130.0   Median :  0.000   Median : 1.600  
##  Mean   : 10.10   Mean   : 86.3   Mean   :  7.556   Mean   : 2.355  
##  3rd Qu.: 13.68   3rd Qu.:130.0   3rd Qu.:  4.000   3rd Qu.: 3.100  
##  Max.   :159.80   Max.   :140.0   Max.   :220.000   Max.   :80.000  
##  NA's   :249                      NA's   :1         NA's   :1       
##    stemheight       livingcanopy     inplotcanopy   
##  Min.   : -1.600   Min.   :  0.00   Min.   : 40.00  
##  1st Qu.:  1.000   1st Qu.: 60.00   1st Qu.:100.00  
##  Median :  1.800   Median :100.00   Median :100.00  
##  Mean   :  4.786   Mean   : 74.25   Mean   : 99.78  
##  3rd Qu.:  4.450   3rd Qu.:100.00   3rd Qu.:100.00  
##  Max.   :134.000   Max.   :100.00   Max.   :100.00  
##                    NA's   :85       NA's   :84

## Continuous variable summaries- Std Devs
sapply(SOAP.cont, sd, na.rm=TRUE)

##             dbh       dbhheight basalcanopydiam   maxcanopydiam 
##       14.927547       58.329815       19.304295        3.803459 
##      stemheight    livingcanopy    inplotcanopy 
##        8.526687       41.145058        2.863307

# matrix scatterplot
plot(SOAP.cont)

#side by side boxplots
boxplot(SOAP.cont, las=2)

Exploratory Analysis for SOAP by Sampling Plots

## by species types by PlotID

xtabs(~plotid+scientificname, data=SOAP.cat)

##           scientificname
## plotid     Abies concolor Amelanchier utahensis Arctostaphylos viscida
##   SOAP1343              0                     0                     12
##   SOAP139               0                     0                     14
##   SOAP143               0                     0                     24
##   SOAP1515              0                     0                     23
##   SOAP1563              0                     0                     14
##   SOAP1611              0                     0                     10
##   SOAP1695              0                     0                     36
##   SOAP187               0                     0                      0
##   SOAP223               0                     0                      8
##   SOAP255               0                     0                     15
##   SOAP283               0                     0                      7
##   SOAP299               0                     0                      3
##   SOAP331               4                     2                      0
##   SOAP43                7                     0                      0
##   SOAP555               0                     0                      6
##   SOAP63                0                     0                      6
##   SOAP95                0                     0                      0
##   SOAP991               0                     0                      1
##           scientificname
## plotid     Calocedrus decurrens Ceanothus cuneatus Ceanothus integerrimus
##   SOAP1343                    1                  5                      0
##   SOAP139                    72                  0                      1
##   SOAP143                    13                  0                     34
##   SOAP1515                    0                  0                      1
##   SOAP1563                    0                  0                      0
##   SOAP1611                    0                  0                      0
##   SOAP1695                    0                  1                      0
##   SOAP187                    81                  0                      0
##   SOAP223                     0                  2                      0
##   SOAP255                    14                  0                      0
##   SOAP283                     3                  0                      1
##   SOAP299                     3                  0                      0
##   SOAP331                    85                  0                      0
##   SOAP43                     63                  0                      0
##   SOAP555                     0                  0                      0
##   SOAP63                     38                  0                      0
##   SOAP95                     61                  0                      5
##   SOAP991                     1                 41                      0
##           scientificname
## plotid     Cercocarpus montanus var. glaber Corylus cornuta
##   SOAP1343                                2               0
##   SOAP139                                 0               0
##   SOAP143                                 0               0
##   SOAP1515                                0               0
##   SOAP1563                                0               0
##   SOAP1611                                0               0
##   SOAP1695                                0               0
##   SOAP187                                 1               0
##   SOAP223                                16               0
##   SOAP255                                 0               0
##   SOAP283                                 0               0
##   SOAP299                                 0               0
##   SOAP331                                 0               0
##   SOAP43                                  0               3
##   SOAP555                                 0               0
##   SOAP63                                  0               0
##   SOAP95                                  0               0
##   SOAP991                                21               0
##           scientificname
## plotid     Pinus lambertiana Pinus ponderosa Prunus emarginata
##   SOAP1343                 0               1                 0
##   SOAP139                  3               6                 0
##   SOAP143                  0              38                 0
##   SOAP1515                 0               1                 0
##   SOAP1563                 0               0                 0
##   SOAP1611                 0               0                 0
##   SOAP1695                 0               4                 0
##   SOAP187                  0               8                 0
##   SOAP223                  0               0                 0
##   SOAP255                  0              51                 0
##   SOAP283                  0               0                 0
##   SOAP299                  0              16                 0
##   SOAP331                  1              11                 0
##   SOAP43                   9               6                 0
##   SOAP555                  0               0                 0
##   SOAP63                   0              53                 0
##   SOAP95                   3              18                 0
##   SOAP991                  0               0                 1
##           scientificname
## plotid     Quercus chrysolepis Quercus garryana Quercus kelloggii
##   SOAP1343                   0                1                 1
##   SOAP139                   19                0                 2
##   SOAP143                    0                0                 1
##   SOAP1515                   1                0                 0
##   SOAP1563                   0                0                 0
##   SOAP1611                   0                0                 0
##   SOAP1695                  50                0                 0
##   SOAP187                    2                0                 2
##   SOAP223                   14                1                 1
##   SOAP255                   15                0                 9
##   SOAP283                    3                0                 0
##   SOAP299                    4                0                 0
##   SOAP331                    6                0                23
##   SOAP43                     0                0                 9
##   SOAP555                    0                0                 0
##   SOAP63                    10                0                15
##   SOAP95                    17                0                 2
##   SOAP991                    0                1                 0
##           scientificname
## plotid     Rhamnus ilicifolia Ribes roezlii
##   SOAP1343                  0             1
##   SOAP139                   0             5
##   SOAP143                   0             4
##   SOAP1515                  0             0
##   SOAP1563                  0             0
##   SOAP1611                  0             0
##   SOAP1695                 12             3
##   SOAP187                   0             1
##   SOAP223                   0             1
##   SOAP255                   0             0
##   SOAP283                   0             0
##   SOAP299                   0             0
##   SOAP331                   0             0
##   SOAP43                    0             3
##   SOAP555                   0             0
##   SOAP63                    0             0
##   SOAP95                    0             0
##   SOAP991                   0             1

## aggregate summary stats 
#  Re-instate PlotID

SOAP.cont$plotid<-insitu.SOAP$plotid

# group by
melted <- melt(SOAP.cont, id.vars=c("plotid"))

SOAP.aggregate<-summarise(group_by(melted, plotid, variable),
          mean=mean(value, na.rm=TRUE), sd=sd(value, na.rm=TRUE), 
          min=min(value, na.rm=TRUE), max=max(value, na.rm=TRUE))
 
SOAP.aggregate

## Source: local data frame [126 x 6]
## Groups: plotid [?]
## 
##      plotid        variable       mean        sd   min   max
##       (chr)          (fctr)      (dbl)     (dbl) (dbl) (dbl)
## 1  SOAP1343             dbh   8.400000  6.437391   1.0  15.4
## 2  SOAP1343       dbhheight 120.000000 33.879582  10.0 130.0
## 3  SOAP1343 basalcanopydiam  29.287500 30.840844   0.0 140.0
## 4  SOAP1343   maxcanopydiam   3.920833  5.736684   0.6  30.0
## 5  SOAP1343      stemheight   2.112500  1.212368   0.7   6.5
## 6  SOAP1343    livingcanopy  71.625000 33.217089   0.0 100.0
## 7  SOAP1343    inplotcanopy 100.000000  0.000000 100.0 100.0
## 8   SOAP139             dbh   7.274510 19.207526   0.4 159.8
## 9   SOAP139       dbhheight  80.983607 59.089648  10.0 130.0
## 10  SOAP139 basalcanopydiam   5.311475 13.882674   0.0  60.0
## ..      ...             ...        ...       ...   ...   ...

## WHOA-- That's a big spicy meatball!
## Break it down by plot again-- this time stratify by parameter

#overside columnn width restriction
options(dplyr.width = Inf)

# dbh
 insitu_dbh <- SOAP.cont %>%
   group_by(plotid) %>%
   summarise(dbh.min = min(dbh, na.rm=TRUE), dbh.max = max(dbh, na.rm=TRUE), dbh.mean=mean(dbh, na.rm=TRUE), dbh.median=median(dbh, na.rm=TRUE),dbh.sd=sd(dbh, na.rm=TRUE))

 insitu_dbh             

## Source: local data frame [18 x 6]
## 
##      plotid dbh.min dbh.max  dbh.mean dbh.median    dbh.sd
##       (chr)   (dbl)   (dbl)     (dbl)      (dbl)     (dbl)
## 1  SOAP1343     1.0    15.4  8.400000       8.60  6.437391
## 2   SOAP139     0.4   159.8  7.274510       1.70 19.207526
## 3   SOAP143     0.4    44.0  6.018889       1.40 11.503276
## 4  SOAP1515     1.0    24.0 11.646154      10.50  6.357753
## 5  SOAP1563     0.0     0.0  0.000000       0.00  0.000000
## 6  SOAP1611     0.0     0.0  0.000000       0.00  0.000000
## 7  SOAP1695     0.0    60.7  7.654667       1.90 12.417656
## 8   SOAP187     1.2   102.9  8.220652       2.75 16.507677
## 9   SOAP223     0.6   107.9 16.823529      10.20 26.028387
## 10  SOAP255     0.4    78.3 10.751724       7.30 13.063271
## 11  SOAP283     0.0    33.5  6.492857       0.00 12.273578
## 12  SOAP299     1.1    64.1 32.685000      33.90 21.016767
## 13  SOAP331     0.4    72.1 10.714729       5.90 12.610971
## 14   SOAP43     0.7    84.0 11.066316       5.80 16.032312
## 15  SOAP555     0.0     0.0  0.000000       0.00  0.000000
## 16   SOAP63     0.2    50.8 14.746364      13.25 12.035816
## 17   SOAP95     0.3    61.3 10.762626       5.50 12.969474
## 18  SOAP991     1.4    10.5  6.420000       8.20  4.490212

## functionalization might be best here...but for now 
 ## JV scripting

 
  # dbhheight
 insitu_dbhheight <- SOAP.cont %>%
   group_by(plotid) %>%
   summarise(dbhheight.min = min(dbhheight, na.rm=TRUE), dbhheight.max = max(dbhheight, na.rm=TRUE), dbhheight.mean=mean(dbhheight, na.rm=TRUE), dbhheight.median=median(dbhheight, na.rm=TRUE),dbhheight.sd=sd(dbhheight, na.rm=TRUE))
 
 insitu_dbhheight

## Source: local data frame [18 x 6]
## 
##      plotid dbhheight.min dbhheight.max dbhheight.mean dbhheight.median
##       (chr)         (dbl)         (dbl)          (dbl)            (dbl)
## 1  SOAP1343          10.0           130      120.00000            130.0
## 2   SOAP139          10.0           130       80.98361            130.0
## 3   SOAP143          10.0           130       59.47368             10.0
## 4  SOAP1515           1.0           130       70.88462            130.0
## 5  SOAP1563           0.0             0        0.00000              0.0
## 6  SOAP1611           0.0             0        0.00000              0.0
## 7  SOAP1695           0.0           130       66.36604             38.5
## 8   SOAP187           2.0           130       69.07368             10.0
## 9   SOAP223          10.0           130       94.18605            130.0
## 10  SOAP255          10.0           140      109.32692            130.0
## 11  SOAP283           0.0           130       28.18571              0.0
## 12  SOAP299          10.0           130      125.38462            130.0
## 13  SOAP331          10.0           130       83.65909            130.0
## 14   SOAP43          10.0           130      103.60000            130.0
## 15  SOAP555           0.0             0        0.00000              0.0
## 16   SOAP63           0.8           130      111.23607            130.0
## 17   SOAP95          10.0           130       85.84906            130.0
## 18  SOAP991          10.0           130      121.04478            130.0
##    dbhheight.sd
##           (dbl)
## 1      33.87958
## 2      59.08965
## 3      59.33022
## 4      65.13944
## 5       0.00000
## 6       0.00000
## 7      61.59260
## 8      60.22498
## 9      55.12896
## 10     45.67221
## 11     55.19133
## 12     23.53394
## 13     58.62428
## 14     49.95998
## 15      0.00000
## 16     43.87607
## 17     58.14379
## 18     31.77260

 # basalcanopydiam
 insitu_basalcanopydiam <- SOAP.cont %>%
   group_by(plotid) %>%
   summarise(basalcanopydiam.min = min(basalcanopydiam, na.rm=TRUE), basalcanopydiam.max = max(basalcanopydiam, na.rm=TRUE), basalcanopydiam.mean=mean(basalcanopydiam, na.rm=TRUE), basalcanopydiam.median=median(basalcanopydiam, na.rm=TRUE),basalcanopydiam.sd=sd(basalcanopydiam, na.rm=TRUE))
 
 insitu_basalcanopydiam

## Source: local data frame [18 x 6]
## 
##      plotid basalcanopydiam.min basalcanopydiam.max basalcanopydiam.mean
##       (chr)               (dbl)               (dbl)                (dbl)
## 1  SOAP1343                 0.0                 140           29.2875000
## 2   SOAP139                 0.0                  60            5.3114754
## 3   SOAP143                 0.0                  40            2.3771930
## 4  SOAP1515                 0.0                 117           23.8461538
## 5  SOAP1563                 2.0                  90           27.2142857
## 6  SOAP1611                 1.1                 220           42.6100000
## 7  SOAP1695                 0.0                  94           10.3000000
## 8   SOAP187                 0.0                  48            1.0315789
## 9   SOAP223                 0.0                  58           11.3255814
## 10  SOAP255                 0.0                  80            5.5000000
## 11  SOAP283                 0.0                  31           14.0000000
## 12  SOAP299                 0.0                 125           15.9230769
## 13  SOAP331                 0.0                  60            0.5174242
## 14   SOAP43                 0.0                  38            0.8320000
## 15  SOAP555                 7.0                 130           90.8333333
## 16   SOAP63                 0.0                  61            2.3360656
## 17   SOAP95                 0.0                  30            0.7809524
## 18  SOAP991                 0.0                 163           34.6567164
##    basalcanopydiam.median basalcanopydiam.sd
##                     (dbl)              (dbl)
## 1                    22.5          30.840844
## 2                     0.0          13.882674
## 3                     0.0           6.012362
## 4                    20.0          21.502916
## 5                    24.0          22.874839
## 6                    18.0          64.730286
## 7                     0.0          17.087255
## 8                     0.0           6.104502
## 9                     6.0          15.085030
## 10                    0.0          15.270379
## 11                   13.0          11.948608
## 12                    0.0          35.305436
## 13                    0.0           5.253389
## 14                    0.0           4.755664
## 15                  110.0          48.803347
## 16                    0.0           8.759987
## 17                    0.0           3.855518
## 18                   24.0          31.471754

 # maxcanopydiam
 insitu_maxcanopydiam <- SOAP.cont %>%
   group_by(plotid) %>%
   summarise(maxcanopydiam.min = min(maxcanopydiam, na.rm=TRUE), maxcanopydiam.max = max(maxcanopydiam, na.rm=TRUE), maxcanopydiam.mean=mean(maxcanopydiam, na.rm=TRUE), maxcanopydiam.median=median(maxcanopydiam, na.rm=TRUE),maxcanopydiam.sd=sd(maxcanopydiam, na.rm=TRUE))
 
 insitu_maxcanopydiam

## Source: local data frame [18 x 6]
## 
##      plotid maxcanopydiam.min maxcanopydiam.max maxcanopydiam.mean
##       (chr)             (dbl)             (dbl)              (dbl)
## 1  SOAP1343               0.6              30.0           3.920833
## 2   SOAP139               0.2              30.0           1.752459
## 3   SOAP143               0.2              11.5           1.240351
## 4  SOAP1515               0.7               5.4           2.334615
## 5  SOAP1563               0.5               3.9           1.885714
## 6  SOAP1611               1.0               5.1           2.240000
## 7  SOAP1695               0.1              60.0           3.513208
## 8   SOAP187               0.0              80.0           2.309474
## 9   SOAP223               0.3              10.1           2.576744
## 10  SOAP255               0.0              12.6           3.031731
## 11  SOAP283               0.0               3.8           1.807143
## 12  SOAP299               0.0              10.1           4.330769
## 13  SOAP331               0.0               8.5           1.934351
## 14   SOAP43               0.0              18.1           2.440000
## 15  SOAP555               1.3               4.0           2.200000
## 16   SOAP63               0.0               6.8           2.604918
## 17   SOAP95               0.0              12.4           2.096226
## 18  SOAP991               0.4               4.7           1.962687
##    maxcanopydiam.median maxcanopydiam.sd
##                   (dbl)            (dbl)
## 1                  2.85        5.7366843
## 2                  0.90        3.2085684
## 3                  0.60        1.8008709
## 4                  2.10        1.0403624
## 5                  1.80        1.0021954
## 6                  1.80        1.3150412
## 7                  2.20        7.1522634
## 8                  0.80        8.3154186
## 9                  1.90        2.1445209
## 10                 2.40        2.3197449
## 11                 1.70        1.2256642
## 12                 4.25        2.9714332
## 13                 1.40        1.6596928
## 14                 1.95        2.4516743
## 15                 1.95        0.9423375
## 16                 2.45        1.6878688
## 17                 0.90        2.2582304
## 18                 2.00        0.9766727

  # stemheight
 insitu_stemheight <- SOAP.cont %>%
   group_by(plotid) %>%
   summarise(stemheight.min = min(stemheight, na.rm=TRUE), stemheight.max = max(stemheight, na.rm=TRUE), stemheight.mean=mean(stemheight, na.rm=TRUE), stemheight.median=median(stemheight, na.rm=TRUE),stemheight.sd=sd(stemheight, na.rm=TRUE))
 
 insitu_stemheight

## Source: local data frame [18 x 6]
## 
##      plotid stemheight.min stemheight.max stemheight.mean
##       (chr)          (dbl)          (dbl)           (dbl)
## 1  SOAP1343            0.7            6.5       2.1125000
## 2   SOAP139            0.5          120.0       4.7459016
## 3   SOAP143            0.5           19.7       2.5307018
## 4  SOAP1515            0.6            3.4       1.7769231
## 5  SOAP1563            0.0            1.9       1.0714286
## 6  SOAP1611            0.6            4.2       1.7200000
## 7  SOAP1695            0.6           17.0       3.0566038
## 8   SOAP187            0.5          134.0       4.8789474
## 9   SOAP223            0.5            8.2       2.4441860
## 10  SOAP255            0.0           44.7       5.1019231
## 11  SOAP283            0.0           13.8       3.3285714
## 12  SOAP299            1.0           33.1      13.8961538
## 13  SOAP331            0.5           39.2       5.0098485
## 14   SOAP43            0.5           51.1       5.0470000
## 15  SOAP555            0.8            1.4       0.9666667
## 16   SOAP63            0.6           33.0       9.2327869
## 17   SOAP95            0.5           28.1       6.0509434
## 18  SOAP991           -1.6            3.8       1.8537313
##    stemheight.median stemheight.sd
##                (dbl)         (dbl)
## 1               2.05     1.2123683
## 2               1.50    13.0803565
## 3               1.00     4.1863057
## 4               1.70     0.5867249
## 5               1.20     0.5426836
## 6               1.30     1.2318008
## 7               2.20     2.8116770
## 8               1.40    14.9537862
## 9               1.60     1.9816182
## 10              2.35     6.7463328
## 11              1.80     4.1795157
## 12             10.30    12.8098550
## 13              2.40     6.4031226
## 14              2.95     7.8090010
## 15              0.90     0.2250926
## 16              5.00     9.0170403
## 17              1.70     8.2020583
## 18              1.90     0.9241422

  # livingcanopy
 insitu_livingcanopy <- SOAP.cont %>%
   group_by(plotid) %>%
   summarise(livingcanopy.min = min(livingcanopy, na.rm=TRUE), livingcanopy.max = max(livingcanopy, na.rm=TRUE), livingcanopy.mean=mean(livingcanopy, na.rm=TRUE), livingcanopy.median=median(livingcanopy, na.rm=TRUE),livingcanopy.sd=sd(livingcanopy, na.rm=TRUE))
  
 insitu_livingcanopy

## Source: local data frame [18 x 6]
## 
##      plotid livingcanopy.min livingcanopy.max livingcanopy.mean
##       (chr)            (int)            (int)             (dbl)
## 1  SOAP1343                0              100          71.62500
## 2   SOAP139                0              100          79.71311
## 3   SOAP143                0              100          95.92105
## 4  SOAP1515               20              100          82.42308
## 5  SOAP1563               NA               NA               NaN
## 6  SOAP1611               NA               NA               NaN
## 7  SOAP1695                0              100          88.57576
## 8   SOAP187                0              100          23.15789
## 9   SOAP223                0              100          79.48837
## 10  SOAP255                0              100          66.89320
## 11  SOAP283               NA               NA               NaN
## 12  SOAP299                0              100          84.61538
## 13  SOAP331                0              100          78.78788
## 14   SOAP43                0              100          81.90000
## 15  SOAP555               NA               NA               NaN
## 16   SOAP63                0              100          62.58197
## 17   SOAP95                0              100          84.90566
## 18  SOAP991                0              100          71.41791
##    livingcanopy.median livingcanopy.sd
##                  (dbl)           (dbl)
## 1                   80        33.21709
## 2                  100        39.42405
## 3                  100        18.47703
## 4                   90        20.29418
## 5                   NA             NaN
## 6                   NA             NaN
## 7                  100        19.37172
## 8                    0        42.40793
## 9                   90        24.98225
## 10                 100        46.67953
## 11                  NA             NaN
## 12                 100        36.79465
## 13                 100        41.03676
## 14                 100        38.57814
## 15                  NA             NaN
## 16                 100        48.22568
## 17                 100        35.96944
## 18                  80        26.73957

   # inplotcanopy
 insitu_inplotcanopy <- SOAP.cont %>%
   group_by(plotid) %>%
   summarise(inplotcanopy.min = min(inplotcanopy, na.rm=TRUE), inplotcanopy.max = max(inplotcanopy, na.rm=TRUE), inplotcanopy.mean=mean(inplotcanopy, na.rm=TRUE), inplotcanopy.median=median(inplotcanopy, na.rm=TRUE),inplotcanopy.sd=sd(inplotcanopy, na.rm=TRUE))
 
 insitu_inplotcanopy 

## Source: local data frame [18 x 6]
## 
##      plotid inplotcanopy.min inplotcanopy.max inplotcanopy.mean
##       (chr)            (int)            (int)             (dbl)
## 1  SOAP1343              100              100         100.00000
## 2   SOAP139               90              100          99.91803
## 3   SOAP143              100              100         100.00000
## 4  SOAP1515              100              100         100.00000
## 5  SOAP1563               NA               NA               NaN
## 6  SOAP1611               NA               NA               NaN
## 7  SOAP1695               65              100          98.03030
## 8   SOAP187              100              100         100.00000
## 9   SOAP223              100              100         100.00000
## 10  SOAP255              100              100         100.00000
## 11  SOAP283               NA               NA               NaN
## 12  SOAP299              100              100         100.00000
## 13  SOAP331              100              100         100.00000
## 14   SOAP43               40              100          98.90000
## 15  SOAP555               NA               NA               NaN
## 16   SOAP63              100              100         100.00000
## 17   SOAP95              100              100         100.00000
## 18  SOAP991              100              100         100.00000
##    inplotcanopy.median inplotcanopy.sd
##                  (dbl)           (dbl)
## 1                  100       0.0000000
## 2                  100       0.9053575
## 3                  100       0.0000000
## 4                  100       0.0000000
## 5                   NA             NaN
## 6                   NA             NaN
## 7                  100       6.7867965
## 8                  100       0.0000000
## 9                  100       0.0000000
## 10                 100       0.0000000
## 11                  NA             NaN
## 12                 100       0.0000000
## 13                 100       0.0000000
## 14                 100       7.7713538
## 15                  NA             NaN
## 16                 100       0.0000000
## 17                 100       0.0000000
## 18                 100       0.0000000

Let’s Extract Some Vegetation Indices for SOAP

## import NDVI
SOAP.ndvi <- raster("../NEONdata/D17-California/SOAP/2013/spectrometer/veg_index/SOAP_NDVI.tif")
SOAP.ndvi

## class       : RasterLayer 
## dimensions  : 1516, 3292, 4990672  (nrow, ncol, ncell)
## resolution  : 1, 1  (x, y)
## extent      : 296906, 300198, 4100038, 4101554  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /Users/R-Files/Documents/data/NEONDI-2016/NEONdata/D17-California/SOAP/2013/spectrometer/veg_index/SOAP_NDVI.tif 
## names       : SOAP_NDVI 
## values      : -0.5165877, 1  (min, max)

##visualize distribution

hist(SOAP.ndvi, main="NDVI for SOAP Site")

## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...):
## 2% of the raster cells were used. 100000 values used.

# visualize raster
plot(SOAP.ndvi, main="NDVI for SOAP Site")

## read in plot centroids
SOAP_plots <- readOGR("../NEONdata/D17-California/SOAP/vector_data/",
                      "SOAP_centroids")

## OGR data source with driver: ESRI Shapefile 
## Source: "../NEONdata/D17-California/SOAP/vector_data/", layer: "SOAP_centroids"
## with 50 features
## It has 10 fields

# Overlay the centroid points and the stem locations on the NDVI plot
plot(SOAP.ndvi,
     main="SOAP Plot Locations \n(N=18)",
     col=gray.colors(100, start=.3, end=.9))

# pch 0 = square
plot(SOAP_plots,
     pch = 0,
     cex = 2,
     col = 4,
     add=TRUE)

# Insitu sampling took place within 40m x 40m square plots, so we use a 20m radius.
# Note that below will return a dataframe containing the max height
# calculated from all pixels in the buffer for each plot (NEON Website)

## we will use these values to populate in-situ data set

NDVI.max <- extract(SOAP.ndvi,
                    SOAP_plots,
                    buffer = 20,
                    fun=(max),
                    sp=TRUE,
                    stringsAsFactors=FALSE)


# NDVI distribution visualization for each plot
# created histograms for all 14 plots of data

cent_ovrList <- extract(SOAP.ndvi,SOAP_plots,buffer = 20)

for (i in 1:14) {
  hist(cent_ovrList[[i]], main=(paste("plot",i)))
  }

Extract descriptive stats from Insitu Data

# import the centroid data and the vegetation structure data

## Call out insitu.SOAP dataframe
## refamilize with PlotIDs within SOAP
unique(insitu.SOAP$plotid)

##  [1] "SOAP43"   "SOAP331"  "SOAP139"  "SOAP1343" "SOAP143"  "SOAP63"  
##  [7] "SOAP1563" "SOAP1695" "SOAP255"  "SOAP1611" "SOAP283"  "SOAP1515"
## [13] "SOAP223"  "SOAP555"  "SOAP299"  "SOAP991"  "SOAP95"   "SOAP187"

# find the max parameter values for each in situ plot
insitu.max <- insitu.SOAP %>%
  group_by(plotid) %>%
  summarise(max.dbh = max(dbh, na.rm=TRUE), max.dbhheight=max(dbhheight, na.rm=TRUE),
            max.basalcanopydiam= max(basalcanopydiam, na.rm=TRUE), 
            max.maxcanopydiam=max(maxcanopydiam), 
            max.stemheight=max(stemheight, na.rm=TRUE),
            max.livingcanopy=max(livingcanopy, na.rm=TRUE),
            max.inplotcanopy= max(inplotcanopy, na.rm=TRUE))

## merge data into data frame, can use approach to create spatial object
## But purposes of capstone will conduct "boring regression analysis"
## merge grouped data (insitu.max) and extracted NDVI (NDVI.max) 
## using plotid to merge

## first must append  NDVI.max file with SOAP

#create data frame for NDVI.max

NDVI.max.df<-NDVI.max@data

# add SOAP prefix to IDs

NDVI.max.df$plot.id<- "SOAP"
NDVI.max.df$plot.id<-(sapply(NDVI.max.df$plot.id, paste0, NDVI.max.df$ID))[,1]

# rename column

colnames(NDVI.max.df)[colnames(NDVI.max.df) == 'SOAP'] <- 'Plot.ID'

SOAP.grand<-merge(NDVI.max.df, insitu.max, 
      by.x = 'plot.id',
      by.y = 'plotid')

# F* Yes

## Reponse variables: dbh, basalcanopydiam, maxcanopydiam, stemheight, livingcanopy, inplotcanopy

## Predictor variable: SOAP_NDVI

## Simple Linear Regression

## Max NDVI vs Max DBH
mod.ndvi<-lm(max.dbh ~ SOAP_NDVI, data=SOAP.grand)
plot(SOAP.grand$max.dbh, SOAP.grand$SOAP_NDVI, xlab="Max NDVI", ylab=" DBH (cm)")
abline(mod.ndvi)

## Max NDVI vs Max Base Canopy Diam
mod2.ndvi<-lm(max.basalcanopydiam~ SOAP_NDVI, data=SOAP.grand)
plot(SOAP.grand$max.basalcanopydiam, SOAP.grand$SOAP_NDVI, 
     xlab="Max NDVI", ylab=" Base Canopy Diameter (cm)")
abline(mod2.ndvi)

## Max NDVI vs Max of the Max Canopy Diam
mod3.ndvi<-lm(max.maxcanopydiam~SOAP_NDVI, data=SOAP.grand)
plot(SOAP.grand$max.maxcanopydiam, SOAP.grand$SOAP_NDVI, 
     xlab="Max NDVI", ylab=" Maximum canopy diameter(m)")
abline(mod3.ndvi)

## Max NDVI vs Max Stem Height
mod4.ndvi<-lm(max.stemheight~SOAP_NDVI, data=SOAP.grand)
plot(SOAP.grand$max.stemheight, SOAP.grand$SOAP_NDVI, 
     xlab="Max NDVI", ylab=" Stem Height (m)")
abline(mod4.ndvi)

create a qualitative landscape assessment index (completely experimental- BE NICE!)

#subset SOAP.grand data set

s<-subset(SOAP.grand, select = c(plot.id, easting, northing, max.dbh, max.dbhheight, max.basalcanopydiam, max.maxcanopydiam, max.stemheight))

## categorize by quartiles, use quantile values from summary and assign score

# dbh
summary(s$max.dbh)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00   41.38   62.70   65.91   88.72  159.80

dbh.cat <- cut(s$max.dbh, breaks = c(-1, 41.38, 62.70, 88.72, 159.8))
dbh.cat<-unclass(dbh.cat)


#basalcanopydiam
summary(s$max.basalcanopydiam)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   30.00   39.50   59.00   77.83   77.00  220.00

basalcanopydiam.cat<- cut(s$max.basalcanopydiam, breaks= c(0, 39.50, 59.0, 77.00, 222.00))
basalcanopydiam.cat<- unclass(basalcanopydiam.cat)

# maxcanopydiam
summary(s$max.maxcanopydiam)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    3.80    5.95   10.10   17.51   15.25   80.00       1

## impute missing value at Plot SOAP331
s[ 8, 7] <- mean(SOAP.grand$max.maxcanopydiam, na.rm=TRUE)

# determine class quantiles scores
summary(s$max.maxcanopydiam)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   3.800   6.375  10.800  17.510  17.660  80.000

max.maxcanopydiam.cat<-cut(s$max.maxcanopydiam, breaks= c(0, 6.375, 10.80, 17.660, 80.00))
basalcanopydiam.cat<-unclass(basalcanopydiam.cat)

s<-data.frame(dbh.cat,basalcanopydiam.cat,basalcanopydiam.cat) 

tree.score<-rowSums(s)

s$tree.score<-tree.score
s$tree.score.perc<- (tree.score/12)

## add on to SOAP.grand dataset

SOAP.grand<-c(SOAP.grand,s)

boxplot(s$tree.score.perc)

## simple linear regression using tree score ()

lm(tree.score.perc~SOAP_NDVI, data=SOAP.grand)

## 
## Call:
## lm(formula = tree.score.perc ~ SOAP_NDVI, data = SOAP.grand)
## 
## Coefficients:
## (Intercept)    SOAP_NDVI  
##      1.4893      -0.9792